In [1]:
import wisps
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import wisps.simulations as wispsim
import astropy.units as u
import splat.empirical as spe
import splat
Purpose: compare various luminosity functions
In [2]:
%matplotlib inline
In [3]:
baraffe=wispsim.make_systems(model_name='baraffe2003', bfraction=0.2)
saumon=wispsim.make_systems(model_name='saumon2008', bfraction=0.2)
sonora=wispsim.make_systems(model_name='marley2019', bfraction=0.2)
phillips=wispsim.make_systems(model_name='phillips2020', bfraction=0.2)
In [4]:
from astropy.io import ascii
In [5]:
klf=pd.read_csv('/users/caganze/research/wisps/data/kirkpatricklf.txt', delimiter=',')
dntb=ascii.read('/users/caganze/research/wisps/data/daniellalf.txt').to_pandas()
In [6]:
def splat_mag_to_spt(mag):
grid=np.arange(17, 39)
rel=spe.typeToMag(grid, 'MKO_J',reference='dupuy2012')[0]
vals= np.polyval(rel, grid)
spt_sorted_idx=np.argsort(vals)
return np.interp(mag, vals[spt_sorted_idx], grid[spt_sorted_idx])
In [7]:
klf['tfm']=np.mean(np.array([klf.t0.values, klf.tf.values]), axis=0)
In [8]:
klf['spt']=klf.tfm.apply(wispsim.splat_teff_to_spt).apply(round)
In [9]:
klf=klf.applymap(float)
In [10]:
dntb['spt']=dntb.M_J.apply(splat_mag_to_spt)
In [11]:
dntb
Out[11]:
In [12]:
klf
Out[12]:
In [13]:
def ryan_lf(J):
logphi=-0.30 + 0.11*(J-14) + 0.15*(J -14)**2.+ 0.015*(J-14)**3-0.00020*(J-14)**4
return (10**logphi)*(10**-3)
def custom_histogram(things, grid, binsize):
n=[]
for g in grid:
n.append(len(things[np.logical_and(g<=things, things< g+binsize)]))
return np.array(n)
In [14]:
from astropy.io import ascii
In [15]:
jgrid=np.arange(10, 18, .5)
teffgrid=np.arange(50, 4000, 150)
In [ ]:
In [16]:
teffs_bar=baraffe['system_teff']
teffs_saumon=saumon['system_teff']
teffs_sonora=sonora['system_teff']
teffs_phil=phillips['system_teff']
normteff_bar = 0.63*(10**-3)/ len(teffs_bar[np.logical_and(teffs_bar>=1650, teffs_bar <=1800)])
normteff_saumon = 0.63*(10**-3)/ len(teffs_saumon[np.logical_and(teffs_saumon>=1650, teffs_saumon <=1800)])
normteff_sonora= 0.63*(10**-3)/ len(teffs_sonora[np.logical_and(teffs_sonora>=1650, teffs_sonora <=1800)])
normteff_phil=0.63*(10**-3)/ len(teffs_phil[np.logical_and(teffs_phil>=1650, teffs_phil <=1800)])
#jmags=wisps.drop_nan(spe.typeToMag(SIMULATED_DIST['spts'][0], '2MASS J')[0])
#jnorm=6.570*(10**-3)/len(jmags[np.logical_and(jmags>=10.25-.25, jmags <=10.25+.25)])
In [17]:
#sem_emp_phi_j=custom_histogram(jmags, jgrid, .5)*jnorm
baraffe_phi_teff=custom_histogram(teffs_bar, teffgrid, 150)*normteff_bar
saumon_phi_teff=custom_histogram(teffs_saumon, teffgrid, 150)*normteff_saumon
sonora_phi_teff=custom_histogram(teffs_sonora, teffgrid, 150)*normteff_sonora
phil_phi_teff=custom_histogram(teffs_phil, teffgrid, 150)*normteff_phil
In [18]:
#count how many things are in the masses of 0.1 and stuff and compare to 0.005 pc to 0.0037 pc^3
#2 things doing this
#teff-> bolometric correction -> simulated luminosity ---> magnitude
#try will best's 2018 relation
#and that we dont't complete samples for some magnitudes
#compare to bochanski's measurement
#look at the scale factors between those two plots you might you want
#you should be plotting the fits that are
#invert relations
#invert
#binary fraction (how many)
#binary fraction random 20%
#binary mass ratio from a distribution from splat (allen et al. from splat)
#the secondary have a magnitude
#the bianay
#mag of the system by combining the flux
#hst magnitude of the secondary by adding the
#adding the
#everything is laid out in burgasser 2007
#educational exercise to find thre number of stars and brown dwarfs
In [19]:
import seaborn as sns
In [20]:
sns.set_palette('viridis')
In [21]:
fig, ax1=plt.subplots(figsize=(8, 5), ncols=1)
#ax.step(jgrid, sem_emp_phi_j, color='#0074D9', label='Simulated')
#ax.step(jgrid, ryan_lf(jgrid), color='#FF4136', label='RyanJr2017')
#ax.errorbar(dntb.M_J, dntb.Density*(10**-3), fmt='o', c='k', label='BG2019')
#ax.set_xlabel('J', fontsize=18)
#ax.set_ylabel(r'LF [pc$^{-3}$ mag$^{-1} $]', fontsize=18)
#ax.minorticks_on()
#ax.legend(fontsize=18)
#ax.set_xlim([10, 16])
#ax.set_ylim([0., 0.005])
ax1.step(teffgrid, baraffe_phi_teff, where='mid')
ax1.step(teffgrid, saumon_phi_teff, where='mid')
ax1.step(teffgrid, sonora_phi_teff, where='mid')
ax1.step(teffgrid, phil_phi_teff, where='mid')
for index, row in klf.iterrows():
if row.lf==0.0:
pass
elif row.lfunc==0.0:
ax1.errorbar(row.tfm, row.lf*(10**-3), yerr=0.0005, color='#111111', fmt='o',lolims=True, ls='none')
else:
ax1.errorbar(row.tfm, row.lf*(10**-3), yerr=row.lfunc*(10**-3), color='#111111', fmt='o')
ax1.set_xlabel(r'T$_{eff}$ [K]', fontsize=18)
ax1.set_ylabel(r'$\frac{dN}{dTeff}$ [ K$^{-1}$ pc$^{-3}$]', fontsize=18)
ax1.minorticks_on()
ax1.legend(fontsize=18, labels=['B03', 'SM03', 'Sonora', 'P20', 'K19',])
ax1.set_xlim([3000, 100.])
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/lfs_teffs.pdf')
In [22]:
js_saumon=wisps.absolute_magnitude_jh(np.sort(saumon['system_spts'].flatten()))[0]
js_bar=wisps.absolute_magnitude_jh(np.sort( baraffe['system_spts'].flatten()))[0]
js_sonora=wisps.absolute_magnitude_jh(np.sort( sonora['system_spts'].flatten()))[0]
js_phil=wisps.absolute_magnitude_jh(np.sort( phillips['system_spts'].flatten()))[0]
In [23]:
jgrid=np.arange(10, 16, .5)
In [24]:
normjs_bar = 1.5*(10**-3)/ len(js_bar[np.logical_and(js_bar>=11.75, js_bar <11.75+0.5)])
normjs_saumon = 1.5*(10**-3)/ len(js_saumon[np.logical_and(js_saumon>=11.75,js_saumon<11.75+0.5)])
normjs_sonora=1.5*(10**-3)/ len(js_sonora[np.logical_and(js_sonora>=11.75,js_sonora<11.75+0.5)])
normjs_phil=1.5*(10**-3)/ len(js_phil[np.logical_and(js_phil>=11.75,js_phil<11.75+0.5)])
baraffe_phi_j=custom_histogram(js_bar, jgrid, 0.5)*normjs_bar
saumon_phi_j=custom_histogram(js_saumon, jgrid,0.5)*normjs_saumon
sonora_phi_j=custom_histogram(js_sonora, jgrid,0.5)*normjs_sonora
phil_phi_j=custom_histogram(js_phil, jgrid,0.5)*normjs_phil
In [25]:
DNLF={"J": dntb.M_J.values, 'lf': dntb.Density.values,
"er":[[1.39, 0.39, 0.29, 0.20, 0.18, 0.15, 0.41, 0.16, 0.15], [1.62, 4.14, 0.31, 0.23, 0.20, 0.18, 1.65, 0.18, 0.18]]}
In [26]:
CRUZ={"J":[10.75, 11.25, 11.75, 12.25, 12.75, 13.25, 13.75, 14.25, 14.75],
"lf": [2.38, 1.66, 1.16, 0.83, 0.50, 0.58, 0.50, 0.66, 0.33],
'er':[0.64, 0.37, 0.31, 0.26, 0.20, 0.22, 0.20, 0.23, 0.17]}
In [27]:
fig, ax1=plt.subplots(figsize=(8, 6), ncols=1)
plt.step(jgrid, ryan_lf(jgrid))
plt.step(jgrid, saumon_phi_j)
plt.step(jgrid, baraffe_phi_j)
plt.step(jgrid, sonora_phi_j)
plt.step(jgrid, phil_phi_j)
plt.errorbar(DNLF['J'], DNLF['lf']*0.001, yerr=np.array(DNLF['er'])*0.001,color='#111111', fmt='o')
plt.errorbar(CRUZ['J'], np.array(CRUZ['lf'])*0.001, yerr=np.array(CRUZ['er'])*0.001,color='#B10DC9', fmt='o')
ax1.set_xlabel(r'$M_J$ [mag]', fontsize=18)
ax1.set_ylabel(r'$\frac{dN}{dM_J}$ [mag $^{-1}$ pc$^{-3}$]', fontsize=18)
ax1.minorticks_on()
ax1.legend(fontsize=18, labels=['R15', 'SM03', 'B03', 'Sonora', 'P20', 'B19',
'C07'])
plt.tight_layout()
plt.savefig(wisps.OUTPUT_FIGURES+'/lfs_js.pdf')
In [32]:
#save the inertpolated luminosity function
LF={}
In [54]:
for modname in ['saumon2008', 'baraffe2003', 'marley2019', 'phillips2020']:
vals=wispsim.make_systems(model_name=modname, bfraction=0.2)
teffs_bar=vals['system_teff']
sys_spts=vals['system_spts']
normteff_bar = 0.63*(10**-3)/ len(teffs_bar[np.logical_and(teffs_bar>=1650, teffs_bar <=1800)])
#normspt
normspt= ((1650-1800)/(24.57-25.77))*0.63*(10**-3)/len(sys_spts[np.logical_and( sys_spts>=24.58, sys_spts <=25.77)])
spt_grid=wisps.splat_teff_to_spt(teffgrid)
phi_spt=custom_histogram(sys_spts, spt_grid, 1)*normspt
LF[modname]={'spt':spt_grid, 'phi':phi_spt}
In [55]:
import pickle
with open(wisps.OUTPUT_FILES+'/lf.pkl', 'wb') as file:
pickle.dump(LF,file)
In [ ]: